#setwd()

library(foreign)
library(splines)
library(survival)

############################################################################
#Read Stata Data File                                                      #    
############################################################################
terrorism <- read.dta("condfrail_terrorism.dta")
summary(terrorism)

############################################################################
#Examine the Data                                                          #    
############################################################################

#Data for Canada
terrorism[24:46, 1:11]

#Data for Mali
terrorism[1220:1242, 1:11]

#Number of Events in Data
numevents <- table(terrorism$attack)
numevents

#Average Number of Events for Each Country
avg.eventcount <- numevents/112
avg.eventcount

#Frequency of Events in Each Strata
num.strata <- table(terrorism$evnum)
num.strata

#Frequency of Events in Each Strata (Collapsed)
num.strata.collap <- table(terrorism$evnumcollap)
num.strata.collap

#Maximum Number of Events for Each Country
maxnum.country <- cbind(terrorism$ccode,terrorism$total.events)
maxnum.country

############################################################################
#Estimate a Conditional Cox Model in Gap Time                              #
############################################################################
cond.gap.collap <- coxph(Surv(istart,ifinish,attack)~trade+fdi+portfolio+
                     economicdevelopment+partnersdevelopment+size+democracy+
                     governmentcapability+conflict+cluster(ccode)+
                     strata(evnumcollap),method=c("efron"),data=terrorism)

###################################
#Obtain the Results from the Model#
###################################
summary(cond.gap.collap)

#Collect the Estimated Beta's
cond.gap.collap.beta <- matrix(cond.gap.collap$coef,ncol=1)

#Collect the Robust S.E's
cond.gap.collap.rse <- matrix(sqrt(diag(cond.gap.collap$var)),ncol=1)

#################################
#Obtain Log-Likelihood for Model#
#################################
cond.gap.collap.loglike2 <- cond.gap.collap$loglik[[2]]
cond.gap.collap.loglike2

#########################################
#Plot Cumulative Hazards for Each Strata#
#########################################
cond.gap.collap.surv <- survfit(cond.gap.collap)
summary(cond.gap.collap.surv)

plot(cond.gap.collap.surv,fun="cumhaz",col=c("black","black","black","black","black","black","black","black","black","black",
     "black"),lty=c(1,2,3,4,5,6),ylab="Cumulative Hazard",xlab="Time (Years)",ylim=c(0,12),lwd=c(2,2,1,1,1,1,1,1,1,1,1))
leg.txt <- c("Event 1","Event 2","Event 3","Event 4","Event 5","Event 6")
leg.txt1 <- c("Event 7","Event 8","Event 9","Event 10","Events 11-23")
legend(3.5,12,leg.txt,col=c("black","black","black","black","black","black"),lty=c(1,2,3,4,5,6),lwd=c(2,2,1,1,1,1))
legend(11,12,leg.txt1,col=c("black","black","black","black","black"),lty=c(1,2,3,4,5),lwd=c(1,1,1,1,1))
box()

############################################################################
#Estimate a Cox model in Elapsed Time with Shared Frailty from Gamma       #
#Distribution using the EM Algorithm                                       #
############################################################################
frailty.gamma.em <- coxph(Surv(estart,efinish,attack)~trade+fdi+portfolio+
                     economicdevelopment+partnersdevelopment+size+democracy+
                     governmentcapability+conflict+
                     frailty(ccode,distribution="gamma",method="em"),
                     method=c("efron"),data=terrorism)

#Estimate the same model without the frailty term to conduct an LR test for theta
nofrailty <- coxph(Surv(estart,efinish,attack)~trade+fdi+portfolio+
                     economicdevelopment+partnersdevelopment+size+democracy+
                     governmentcapability+conflict,method=c("efron"),data=terrorism)

###################################
#Obtain the Results from the Model#
###################################
summary(frailty.gamma.em)

#Collect the Estimated Beta's
frailty.gamma.em.beta <- matrix(frailty.gamma.em$coef,ncol=1)

#Collect the S.E's
frailty.gamma.em.se <- matrix(sqrt(diag(frailty.gamma.em$var)),ncol=1)

#################################
#Obtain Log-Likelihood for Model#
#################################
frailty.gamma.em.loglike2 <- frailty.gamma.em$loglik[[2]]
frailty.gamma.em.loglike2

################################################
#Obtain Estimated Variance of the Random Effect#
################################################
frailty.gamma.em.theta <- frailty.gamma.em$history[[1]]$theta
frailty.gamma.em.theta

###################
#LR test for theta#
###################
#Obtain I-likelihood
frailty.gamma.em.ill <- frailty.gamma.em$history[[1]]$c.loglik 
frailty.gamma.em.ill

#Obtain likelihood of the no frailty model
nofrailty.loglike2 <- nofrailty$loglik[[2]] #Base likelihood for LLR test
nofrailty.loglike2

#Calculate LR for theta, report this test for whether theta is significant
frailty.gamma.em.lr <- 2*(4245.22-4127.47)
frailty.gamma.em.lr

#Calculate p-value, the test has 1 df
frailty.gamma.em.lrtest <- 1-pchisq(frailty.gamma.em.lr,df=1)
frailty.gamma.em.lrtest

#################################################
#Calculate and Plot Profile Likelihood for theta#
#################################################
pllik.frailty.gamma.em <- double(100)
for(i in 1:100){
  a <- .01*i
  temp2 <- coxph(Surv(estart,efinish,attack)~trade+fdi+portfolio+
                     economicdevelopment+partnersdevelopment+size+democracy+
                     governmentcapability+conflict+
                     frailty(ccode,theta=a),method=c("efron"),data=terrorism)
  pllik.frailty.gamma.em[i] <- temp2$history[[1]]$c.loglik
}

a <- seq(.01,1,length=100)
temp <- frailty.gamma.em$history[[1]]$c.loglik - 3.84/2

#Calculate confidence interval using temp
frailty.gamma.em.ci <- cbind(a,pllik.frailty.gamma.em)

#Plot Profile Likelihood
plot(a, pllik.frailty.gamma.em, type='l', ylab="Log (Partial-Likelihood)",ylim=c(-4200,-4120),xlab=expression(paste(theta)))
abline(h=temp, lty=2)
points(frailty.gamma.em.theta,frailty.gamma.em.ill)
text(x=.56,y=-4125,adj=c(.25,.25),c(expression(paste(hat(theta)))))
arrows(.55,-4125,.525,-4126.5,length=0.025)

########################################
#Obtain Country Level Frailty Estimates#
########################################
frailties.gamma.em <- frailty.gamma.em$frail
summary(frailties.gamma.em)

#Obtain Country Labels
country.abb <- c("USA","CAN","HAI","DOM","JAM","TRI","MEX","GUA","HON","SAL","NIC","COS","PAN","COL","VEN","GUY",
                 "ECU","PER","BRA","BOL","PAR","CHL","ARG","URU","UKG","IRE","NTH","FRN","SWZ","SPN","POR","GMY",
                 "POL","AUS","HUN","ITA","ALB","SLV","GRC","CYP","BUL","ROM","RUS","EST","LAT","LIT","UKR","BLR",
                 "ARM","SWD","NOR","DEN","ICE","MLI","SEN","MAA","NIR","CDI","BFO","SIE","GHA","TOG","CAO","NIG",
                 "GAB","CEN","CHA","CON","KEN","TAZ","BUI","RWA","ETH","ANG","MZM","ZAM","ZIM","MAW","SAF","NAM",
                 "LES","BOT","SWA","MOR","ALG","TUN","IRN","TUR","EGY","SYR","JOR","ISR","SAU","KUW","BAH","CHN",
                 "ROK","JPN","IND","PAK","SRI","NEP","THI","LAO","MAL","SIN","PHI","INS","AUL","PNG","NEW","FIJ")

frailty.gamma.em.estimates.country <- cbind(country.abb, frailties.gamma.em)

########################################
#Plot Frailty Estimate for Each Country#
########################################
plot(frailties.gamma.em,type='n',ylab="Frailty Estimates",xlab="Country Number",ylim=c(-2.2,2.2),xlim=c(0,120))
text(frailties.gamma.em,country.abb)
abline(h=0,col="black",lty=2,lwd=1)
box()

#########################
#Plot Cumulative Hazards#
#########################
frailty.gamma.em.surv <- survfit(frailty.gamma.em,conf.type=c("none"))
summary(frailty.gamma.em.surv)

plot(frailty.gamma.em.surv,fun="cumhaz",col="black",lwd=2, ylab="Cumulative Hazard",xlab="Time (Years)",xlim=c(0,23))
box()

############################################################################
#Estimate a Conditional Frailty Model in Gap Time with Shared Frailty from #
#Gamma Distribution using the EM Algorithm                                 #
############################################################################
condfrailty.gamma.em <- coxph(Surv(istart,ifinish,attack)~trade+fdi+portfolio+
                              economicdevelopment+partnersdevelopment+size+democracy+
                              governmentcapability+conflict+strata(evnumcollap)+
                              frailty(ccode,distribution="gamma",method="em"),
                              method=c("efron"),data=terrorism)

#Estimate the same model without the frailty term to conduct an LR test for theta
condnofrailty <- coxph(Surv(istart,ifinish,attack)~trade+fdi+portfolio+
                              economicdevelopment+partnersdevelopment+size+democracy+
                              governmentcapability+conflict+strata(evnumcollap),
                              method=c("efron"),data=terrorism)

###################################
#Obtain the Results from the Model#
###################################
summary(condfrailty.gamma.em)

#Collect the Estimated Beta's
condfrailty.gamma.em.beta <- matrix(condfrailty.gamma.em$coef,ncol=1)

#Collect the S.E's
condfrailty.gamma.em.se <- matrix(sqrt(diag(condfrailty.gamma.em$var)),ncol=1)

#################################
#Obtain Log-Likelihood for Model#
#################################
condfrailty.gamma.em.loglike2 <- condfrailty.gamma.em$loglik[[2]]
condfrailty.gamma.em.loglike2

################################################
#Obtain Estimated Variance of the Random Effect#
################################################
condfrailty.gamma.em.theta <- condfrailty.gamma.em$history[[1]]$theta
condfrailty.gamma.em.theta

###################
#LR test for theta#
###################
#Obtain I-likelihood
condfrailty.gamma.em.ill <- condfrailty.gamma.em$history[[1]]$c.loglik 
condfrailty.gamma.em.ill

#Obtain likelihood of the no frailty model
condnofrailty.loglike2 <- condnofrailty$loglik[[2]] #Base likelihood for LLR test
condnofrailty.loglike2

#Calculate LR for theta, report this test for whether theta is significant
condfrailty.gamma.em.lr <- 2*(4296.36-4295.09)
condfrailty.gamma.em.lr

#Calculate p-value, the test has 1 df 
condfrailty.gamma.em.lrtest <- 1-pchisq(condfrailty.gamma.em.lr,df=1)
condfrailty.gamma.em.lrtest

#################################################
#Calculate and Plot Profile Likelihood for theta#
#################################################
pllik.condfrailty.gamma.em <- double(100)
for(i in 1:100){
  a <- .01*i
  temp2 <- coxph(Surv(istart,ifinish,attack)~trade+fdi+portfolio+
                              economicdevelopment+partnersdevelopment+size+democracy+
                              governmentcapability+conflict+strata(evnumcollap)+
                              frailty(ccode,theta=a),
                              method=c("efron"),data=terrorism)
  pllik.condfrailty.gamma.em[i] <- temp2$history[[1]]$c.loglik
}

a <- seq(.01,1,length=100)
temp <- condfrailty.gamma.em$history[[1]]$c.loglik - 3.84/2

#Calculate confidence interval using temp
condfrailty.gamma.em.ci <- cbind(a,pllik.condfrailty.gamma.em)

#Plot Profile Likelihood
plot(a, pllik.condfrailty.gamma.em, type='l', ylab="Log(Partial-Likelihood)",xlab=expression(paste(theta)),ylim=c(-4330,-4280))
abline(h=temp, lty=2)
points(condfrailty.gamma.em.theta, condfrailty.gamma.em.ill)
text(x=.08,y=-4292.5,adj=c(.25,.25),c(expression(paste(hat(theta)))))
arrows(.07,-4293,.045,-4294.5,length=0.025)

########################################
#Obtain Country Level Frailty Estimates#
########################################
condfrailties.gamma.em <- condfrailty.gamma.em$frail
summary(condfrailties.gamma.em)

condfrailty.gamma.em.estimates.country <- cbind(country.abb, condfrailties.gamma.em)

########################################
#Plot Frailty Estimate for Each Country#
########################################
plot(condfrailties.gamma.em,type='n',ylab="Frailty Estimates",xlab="Country Number",ylim=c(-.25,.25),xlim=c(0,120))
text(condfrailties.gamma.em,country.abb)
abline(h=0,col="black",lty=2,lwd=2)
box()

#########################################
#Plot Cumulative Hazards for Each Strata#
#########################################
condfrailty.gamma.em.surv <- survfit(condfrailty.gamma.em)
summary(condfrailty.gamma.em.surv)

plot(condfrailty.gamma.em.surv,fun="cumhaz",col=c("black","black","black","black","black","black","black","black",
     "black","black","black"),lty=c(1,2,3,4,5,6),ylab="Cumulative Hazard",xlab="Time (Years)",ylim=c(0,12),lwd=c(2,2,1,1,1,1,1,1,1,1,1))
leg.txt <- c("Event 1","Event 2","Event 3","Event 4","Event 5","Event 6")
leg.txt1 <- c("Event 7","Event 8","Event 9","Event 10","Events 11-23")
legend(2.5,12.0,leg.txt,col=c("black","black","black","black","black","black"),lty=c(1,2,3,4,5,6),lwd=c(2,2,1,1,1,1))
legend(10,12.0,leg.txt1,col=c("black","black","black","black","black"),lty=c(1,2,3,4,5),lwd=c(1,1,1,1,1))
box()

#######################################################################################
#Collect the Results the Conditional Gap Time, Frailty, and Conditional Frailty Models#                                                       
#######################################################################################

results.beta <- cbind(cond.gap.collap.beta,frailty.gamma.em.beta,condfrailty.gamma.em.beta)
print.default(results.beta,2)

results.se <- cbind(cond.gap.collap.rse,frailty.gamma.em.se,condfrailty.gamma.em.se)
print.default(results.se,2)

results.loglike <- cbind(cond.gap.collap.loglike2,frailty.gamma.em.loglike2,condfrailty.gamma.em.loglike2)
print.default(results.loglike,6)

results.Ilike <- cbind(frailty.gamma.em.ill,condfrailty.gamma.em.ill)
print.default(results.Ilike,6)

results.theta <- cbind(frailty.gamma.em.theta,condfrailty.gamma.em.theta)
print.default(results.theta,2)

results.theta.lr <- cbind(frailty.gamma.em.lr,condfrailty.gamma.em.lr)
print.default(results.theta.lr,2)

results.theta.lrtest <- cbind(frailty.gamma.em.lrtest,condfrailty.gamma.em.lrtest)
print.default(results.theta.lrtest,2)
